The Stochastic Spectral Dynamics of Bending and Tumbling 
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Traditional models of wormlike chains in shear flows at finite temperature approximate the equa- 
tion of motion via finite difference discretization (bead and rod models) . We introduce here a new 
method based on a spectral representation in terms of the natural eigenfunctions. This formula- 
tion separates tumbling and bending dynamics, clearly showing their interrelation, naturally orders 
the bending dynamics according to the characteristic decay rate of its modes, and displays coupling 
among bending modes in a general flow. This hierarchy naturally yields a low dimensional stochastic 
dynamical system which recovers and extends previous numerical results and which leads to a fast 
and efficient numerical method for studying the stochastic nonlinear dynamics of semiflexible poly- 
mers in general flows. This formulation will be useful for studying other physical systems described 
by constrained stochastic partial differential equations. 
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In fields such as complex fluids and single-molecule 
biophysics, one is often interested in dynamic and 
statistical properties of conformations of semiflexible 
biopolymers, those for which resistance to bending is 
comparable to or larger than thermal forcing. These 
polymers may be considered as "wrinkled" rods, in 
contrast with flexible polymers whose coiled config- 
urations are modeled as diffusive random walks in 
space Q. Examples of the two extremes are mi- 
crotubules, responsible (among other roles) for sep- 
arating chromosomes during cell division, and long 
strands of DNA, e.g., those used in the flow experi- 
ments of Chu and coworkers 0] 

Numerical models of semiflexible polymers often rely 
on the bead-spring or bead-rod H, 0, discretiza- 
tions of the wormlike chain model. Such discretiza- 
tion schemes were originally developed for model- 
ing fully flexible polymers as hundreds of statis- 
tically independent extensible or inextensible links 
H, S S @ ■ Modeling becomes progressively more 
complicated as the segments are treated as rods, 
because the inextensibility constraint must be en- 
forced by tensions which are determined by an al- 
gebraic (involving no time derivatives) constraint 
which couples the motion of all the parts of the 
chain. Introducing inextensibility presents two addi- 
tional challenges: multiplicative coupling of the ten- 



sion and the configuration of the polymer [ToL Hl| , 
which leads to multiplicative noise in the polymer 
equation of motion [9j ; and a subtle correlation be- 
tween the configuration of the links which comes 
from projecting out the momenta degrees of free- 
dom of a system which is constrained to live on a 
non-Euclidian subspace 0, 0, 0, Q . Both chal- 
lenges can be met by using appropriate simulation 
algorithms Q — e.g., a mid-step algorithm |rjRlJl3|. 
or an appropriate predictor-corrector algorithm Isl- 
and by introducing appropriate "metric" forces into 
the equations of motion 0, 0, F3 and computing 
them with efficient procedures Introducing a 
bending energy in the model yields further compli- 
cations because inextensibility couples the longitu- 
dinal and transversal chain dynamics. The relax- 
ation time of a bending mode of wavenumber j is 
tj = (±L 4 /(kjK) 15], whereas the longest relax- 
ation time of a rodlike chain is T rot = (±L 3 /(72kBT). 
Here C_l ~ 47t/k/ ln(L/p) is the perpendicular friction 
coefficient, L is the chain length, p is the chain di- 
ameter, [i is the liquid viscosity, kj ~ (2j + l)n is 
the dimensionless relaxation rate of the j-th bend- 
ing eigenmode, k is the chain bending stiffness, ks 
is Boltzmann's constant, and T is temperature. In 
simulations, the fast shape fluctuations of the chain 
must be resolved to sample correctly the Boltzmann 
distribution; therefore, the simulation timestep must 
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be roughly A t » 0.1C±L 4 /((2iV + 1) 4 tt 4 k), where iV 
is the shortest wavenumber that can be captured by 
the discretized model and coincides with the num- 
ber of rods (minus one) in a bead-rod model. Sim- 
ulating the terminal relaxation of a chain while re- 
solving the fast modes requires w r rot /At timesteps, 
i.e, w 10(27V+ l) 4 7r%/(72L), where L^ = n/k B T 
is the persistence length of the chain [ljj. For the 
calculation of linear viscoelastic properties, such a 
hurdle was overcome by simulating the short-time 
behavior of semiflexiblc molecules with finely dis- 
cretized chains, and collapsing the results with those 
obtained by simulating the long-time behavior with 
coarsely-discretized chains |2l Such strategy, 

however, does not work when applied to nonlinear 
viscoelastic properties and general flows and defor- 
mations, because semiflexible chains buckle when 
subjected to strong velocity gradients thus, 
in nonlinear simulations one has to retain the fine 
discretization (in order to resolve curvature during 
buckling) while also simulating long timescales (in 
order to capture the steady rheological behavior). 
Therefore, a new, effective method is needed for 
studying flows of stiff chains at finite temperature. 

As in earlier studies, the starting point is the equa- 
tion of motion for a continuous curve of position 
r(s, t) as a function of arclength s and time t, of 
bending stiffness k, with inextensibility enforced by 
a tension A, moving in a viscous liquid at tempera- 
ture T and velocity U w r ■ VU = jr fillip: 

d t v - -yr = V (9 s F c i astic + <9 s F tc „ s ion + £) (1) 

where (£(s, t)£(s', t')) = 2k B TV~ 1 6(t - t')S{s - s'), 
Fciastic = -nd^r, Ptension = A<9 s r, and the ten- 
sion A is determined by satisfying the inextensi- 
bility constraint \d s r\ = 1. The mobility tensor 
V = ciM + c || tt enforces the anisotropy associ- 
ated with slender bodies in Stokes flow (cj_ = 1/Ci 
and cy = 2c± for a slender cylinder). Unlike in pre- 
vious studies, we work with the Eq. ^ rather than 
discretizing it into a bead-rod formulation. We pose 
the dynamic in two dimensions, separating bending 
from tumbling without constraining the position or 
orientation of the rod whose configuration is the base 
state of our perturbation. Differentiating Eq. ^ and 
recasting the results in terms of tangent and normal 
vectors t = d s r = (cos d, sin??), n = (— sin?9, cos 
yields 

$ t n - 7t = d s V (<9 s F e i ast ic + <9 s F tc n S ion + £) ■ 

Functionally differentiating the elastic energy E = 
(k/2) J L ds(c? s i?) 2 yields the elastic force; whereas 



the orientation i9 does not appear in the energy, the 
curvature d s "d is suppressed; this introduces a natu- 
ral perturbation parameter e 2 = k B TL/ n = L/L p . 
Neglecting terms nonlinear in the curvature i) s or its 
derivatives, 

d t « n-yt(i?) -c ± nd^ + . .. (2) 

suggests representing $ as a combination of eigen- 
functions of the leading order operator: L 4 <9 4 W a = 
fc 4 W Q . These biharmonic eigenfunctions gener- 
alize many properties of Fourier modes and, for 
high mode number, asymptote to the Fourier ba- 
sis 0,0,|2(j. However, they satisfy the boundary 
conditions (no traction and no torque) imposed on 
the ends of an elastic rod immersed in a viscous fluid: 
d s W a = d 2 W a = 0. The only allowed 0-mode is 
independent of s. In this way, the biharmonic eigen- 
function basis naturally yields the representation 

d(t, s) = o (t)W°(s) + $ a (t)W a {s) =$ + $ a W a 

where $o(i) is the rod or tumbling term and the 
"d a (t) are a dynamically ordered sequence of bending 
terms (summation implied). Adjoint to the set of 
shape functions is a set of weight eigenfunctions W a 
of conjugate boundary conditions W a — d^W a = 0. 
The {W a , W a } are biorthogonal / dsW a W fj = L5@, 
a crucial observation for equipartition of energy, 

which yields immediately the equilibrium variance of 
any mode (i? 2 ) = (e/k a ) 2 . This representation of- 
fers a natural way of capturing the long wavelength 
bending modes with accuracy while suppressing the 
short-wavelength modes (which are uninportant in 
hydrodynamics) by truncating at a small a. The 
adjoint functions form a natural basis for the ten- 
sion, which obeys A = at edges: A = A a W a 

The dynamics of the modes and the constraint can 
be derived in simple shear (7 = 7xy T ) by the 
Galerkin method by multiplying Eq.^by the appro- 
priate weight functions (W a for the modes, W Q for 
the tensions), integrating along the chain contour, 
taking the limit for small s, and dropping terms 
0(e 2 ) from the A equation and 0(e 3 ) from the $ 
equation, 

= -c\\kpA^ + + [r)^]p (4) 

tf = -7 (sin 2 tf + cos(2tf )[tf 2 ]o) + + [77% 

d = [n-yn]p - -jj- hp-ftp + t^S^L + [t, x ] p 
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FIG. 1: Theoretical curve (^/ k B T L / {nk 2 ) , Eq. EJ and 
computed values of sj (-da), first 9 bending modes. 



where summation over a, 7 is implied, rj^'^ — 
{t, n} • d s £, projection onto the (3 mode is indicated 
by [• • - ]/3, and 5 is an overlap integral of the ba- 
sis functions which depends on the mobility con- 
stants {cikcll}. For isotropic drag, for example, 

S« 7 oc / dsW (kt,WPW a - fc 2 W a W^ . 

Stability - Shear-induced buckling, in which shear 
and elasticity compete to determine the stability of 
a non-Brownian rod at fixed has been studied an- 
alytically, experimentally |2l|, and numerically |2^ . 
The numerical results show that in simple shear 
the instability occurs when (7»L 4 sin2$o)/(c||K) = 
— 153.2. We derive analytically this criterion by ex- 
panding the shear terms 

7 j_(0) = n 7 n « 7 ± (i?o) + $ a W a i ± (d ) (5) 

(+0($ 2 )) — in simple shear 7_l($o) oc — sin 2 $ - Ex- 
panding the equation of motion and dropping •&- and 
T— dependent terms from the tension, the linear in- 
stability criterion is 



= 



sin 2Qn 



C_1_K, 4 

~JT k P 



tip + (6) 



The overlap integral S^ is well- approximated by 

its diagonal contribution E^ Q 8p. Using the large- 
en Fourier approximation of the basis functions 
yields the instability criterion %L 4 sixi2'3o/(cuk) = 
-1627r 4 /(76 + 37r 2 ) = -149.4, with an error of 2.6%. 
Calculating the diagonal overlap integral reduces the 
error to .2% El. 



FIG. 2: Stochastic spectral dynamics of bending and 
tumbling. The curve which leaves the origin traces #o(£)- 
The remaining curves describe the bending modes $ a (t). 
Each successive mode is smaller in autocorrelation, as 
expected from the equipartion (Eq. [SJ . Here 7 = 8.5 ■ 
10 3 kL 4 /(cx), L — L p , and the x-axis is tj. 



Finite temperature - In shear flow, the zero- 
temperature dynamics asymptote to a straight rod 
aligned in the velocity direction, because the equa- 
tion of change of i?o reduces to #0 = —7 sin 2 ^ — ^ 
as — > jir. In the presence of stochastic forces, 
the dynamics is richer because the $0 = solution 
is nonlinearly unstable; thus, the stochastic forces 
episodically drive $0 to initiate tumbling events, as- 
sisted by shear-induced buckling en route. To study 
the fully-nonlinear stochastic dynamics we use a 
midstep algorithm, consistent with variable diffusiv- 
ity and Stratonovich-interpreted noise |24l |25j . In 
the absence of shear, the algorithm yields the correct 
Boltzmann distribution of bending modes (Fig. EJ. 
Illustrative dynamics in the presence of shear are 
shown in Figs. and 02 

Extensions and connections - One motivation for 
this work is better understanding of complex fluids. 
In dilute suspensions, the contribution of molecules 
to the stress tensor can be calculated via the contin- 
uous generalization for wormlikc chains j2|| of the 
Kirkwood formula |27|. Quantitative comparison 
between the efficient spectral numerics and calcula- 
tions of the stress and molecular conformation (e.g., 
gyration tensor) by established but inefficient large- 
scale simulations of bead-rod chains is in progress 
|ri| . Albeit restricted to two dimensions, the fully 
nonlinear stochastic method is completely general 
with respect to time-dependent flows and kinematics 
(shear, extensional, mixed); thus, it will be useful for 
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FIG. 3: The same dynamics, represented in real space 
and imaged at several representative values of t-y. The 
complex motion is captured by only 6 degrees of freedom. 



comparison with "single-bead rheology" studies in 
actin gels [3] . An extension of the spectral method 
to include the full three-dimensional representation 
of semiflexible chains in general flows is underway. 

An additional topic for future work is the application 
of stochastic mode-elimination 29, 30, 31, 32], inte- 
grating out the fastest bending modes. Mode elimi- 
nation rcnormalizes the noise statistics, and can be 
used to yield a low-dimensional partial differential 
equation for the probability density. Because the 
properties of interest for biopolymers (e.g., configu- 
rational changes connected with macroscopic rheol- 
ogy) evolve on slow timescales, stochastic mode elim- 
ination is an excellent candidate for mathematically 
sound yet experimentally relevant analysis. Such a 
foray would be quite daunting without simple, low 
dimensional stochastic dynamical systems as in Eq.^ 
for which the mapping from Langevin descriptions of 
trajectories to Fokker-Planck descriptions of proba- 
bilities is well-developed. 
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